How many distinct group of population are there in my dataset?

Typical clustering methods (except for some methods like hierarchical clustering) require users to specify number of distinct groups (k) to cluster from your data. In an unsupervised setting, we do not know the exact number and therefore we will run clustering for various number of k.

# We will not run for various `k` to save time. Instead, we will load pre-computed results for `k` between 3 to 8

# This is an easy way to run `scClust` for k = (3,4,5,6,7,8).
# all_k = 3:8
# simlr_results = sapply(as.character(all_k), function(k) {
#   scClust(mat, as.numeric(k), similarity = "pearson", method = "simlr", seed = 1, cores.ratio = 0, geneFilter = 0)
# }, USE.NAMES = TRUE, simplify = FALSE)

# For demonstration purpose, we will run k = 6 (which is actually the number of cell types in our dataset)
simlr_result_k6 = scClust(mat, 6, similarity = "pearson", method = "simlr", seed = 1, cores.ratio = 0, geneFilter = 0)

## Computing the multiple Kernels.
## Performing network diffiusion.
## Iteration:  1 
## Iteration:  2 
## Iteration:  3 
## Iteration:  4 
## Iteration:  5 
## Iteration:  6 
## Iteration:  7 
## Iteration:  8 
## Iteration:  9 
## Iteration:  10 
## Performing t-SNE.
## Epoch: Iteration # 100  error is:  0.1527458 
## Epoch: Iteration # 200  error is:  0.1218641 
## Epoch: Iteration # 300  error is:  0.1112158 
## Epoch: Iteration # 400  error is:  0.1059176 
## Epoch: Iteration # 500  error is:  0.1026344 
## Epoch: Iteration # 600  error is:  0.1003277 
## Epoch: Iteration # 700  error is:  0.0985936 
## Epoch: Iteration # 800  error is:  0.09722343 
## Epoch: Iteration # 900  error is:  0.09610318 
## Epoch: Iteration # 1000  error is:  0.09516611 
## Performing Kmeans.
## Performing t-SNE.
## Epoch: Iteration # 100  error is:  11.56649 
## Epoch: Iteration # 200  error is:  0.2845774 
## Epoch: Iteration # 300  error is:  0.2560054 
## Epoch: Iteration # 400  error is:  0.2410104 
## Epoch: Iteration # 500  error is:  0.2361371 
## Epoch: Iteration # 600  error is:  0.2331395 
## Epoch: Iteration # 700  error is:  0.231009 
## Epoch: Iteration # 800  error is:  0.229431 
## Epoch: Iteration # 900  error is:  0.2281798 
## Epoch: Iteration # 1000  error is:  0.2271669

# load("~/Dropbox (Sydney Uni)/workdir/singleCell/workshop/data/simlr_result_k6.RData")

load(paste(getwd(), "/data/simlr.results.RData", sep = ""))

How do we select optimal k?

If our k is correct, we expect that the each clusters are closely packed together and the distance between the clusters, within-cluster sum of squares (WSS) are expected to be small. Thus, we will select the k with a small total WSS (compact clusters). This is called the “Elbow” method.

# Find total WSS from all cluster outputs
all_wss = sapply(simlr_results, function(result) {
  sum(result$y$withinss)
}, USE.NAMES = TRUE, simplify = TRUE)

plt.dat = data.frame(
  k = names(all_wss),
  total_wss = all_wss
)

ggplot(plt.dat, aes(x = as.numeric(as.character(k)), y = total_wss)) +
  geom_point(size = 3) +
  stat_smooth(method = loess, col = "red", method.args = list(degree = 1)) +
  ylab("Total WSS") +
  ggtitle("Compare Total WSS for each 'k'")

As shown in this plot, the graph begin to plateau from k = 5. We can estimate that k is around 5 or 6. We can further investigate using silhouette scores or other metrics but using our t-SNE plot (and with a little cheating) we will estimate k = 6.

Quick quiz

1. We can determine k with compactness of the clusters. What other measure can we use to determine k?
2. How would you determine optimal number of clusters for hierarchical clustering?

For the purpose of demonstration, we would like to highlight the effect of similarity metric to your cluster output.

# To run scClust with euclidean distance, uncommnet the following lines.
# simlr_result_eucl_k6 = scClust(mat, 6, similarity = "euclidean", method = "simlr", seed = 1, cores.ratio = 0, geneFilter = 0)

# for convenience, we will load our pre-computed result
load(paste(getwd(), "/data/simlr_result_eucl_k6.RData", sep = ""))

# create tsne object
set.seed(123)
tsne_result = Rtsne(t(mat), check_duplicates = FALSE)

#################################################
tmp_lab = as.numeric(factor(lab))
pear_cluster = mapvalues(
  simlr_result_k6$y$cluster,
  from = c(1,2,3,4,5,6),
  to = c(2,3,1,4,6,5)
)
eucl_cluster = mapvalues(
  simlr_result_eucl_k6$y$cluster,
  from = c(1,2,3,4,5,6),
  to = c(6,2,4,3,1,5)
)
#################################################

plt.dat = data.frame(
  tsne1 = rep(tsne_result$Y[,1], 3),
  tsne2 = rep(tsne_result$Y[,2], 3),
  cluster = factor(c(tmp_lab, pear_cluster, eucl_cluster)),
  label = factor(c(rep("Truth", length(lab)), rep("Pearson", length(lab)), rep("Euclidean", length(lab))), levels = c("Truth", "Pearson", "Euclidean"))
)

ggplot(plt.dat, aes(x = tsne1, y = tsne2, colour = cluster, colors)) +
  geom_point(size = 2) +
  ggtitle("t-SNE plot") +
  facet_grid(cols=vars(label)) + 
  theme(legend.position = "none")

Quick quiz

1. From your observation of the above t-SNE plot, did our clustering method (with Pearson correlation) group cells well?

[Optional] We can evaluate the clustering performance using our ground truth.

# ARI
ari = c(mclust::adjustedRandIndex(lab, simlr_result_eucl_k6$y$cluster), mclust::adjustedRandIndex(lab, simlr_result_k6$y$cluster))

# NMI
nmi = c(igraph::compare(as.numeric(factor(lab)), simlr_result_eucl_k6$y$cluster, method = "nmi"), igraph::compare(as.numeric(factor(lab)), simlr_result_k6$y$cluster, method = "nmi"))

plt.dat = data.frame(
  dist = rep(c("Euclidean", "Pearson"), 2),
  value = c(ari, nmi),
  eval = rep(c("ARI", "NMI"), each = 2)
)

ggplot(plt.dat, aes(x = dist, y = value, fill = dist)) + 
  geom_bar(stat="identity") + 
  facet_grid(col = vars(eval)) +
  labs(x = "Similarity metrics", y = "Evalution score", title = "Affect of similarity metrics in scRNA-seq data")

2.1 What are the marker genes that distinguish the different cell types?

Here we provide a function that allow one to find differentially expressed gene between a cluster and the remaining clusters. The input of this function is the cluster id. The output is a list of gene and its associated p - value.

# This function is specific for scClust ouptut object that ran with "SIMLR"
findmarker <- function(mat, scClust_obj, cluster_id){

  #group 1 is the cluster of interest
  group1 = which(scClust_obj$y$cluster == cluster_id)  
  group1 = colnames(mat)[group1]  #find out which cell belongs to group 1
  
  group_label <- ifelse(colnames(mat) %in% group1, "group1", "group2")  #if the cell is in group 1 ,create a label called group 1, else create a label called group 2
  
  #create a dataframe , with cell name and the group that this cell belongs to 
  group_label <- data.frame(group_label, colnames(mat))
  rownames(group_label) <- group_label[,2]
  colnames(group_label) <- c("group","wellKey")
  #MAST requires a wellkey (which essentially is the name of the cell)
  
  
  #create a dataframe, with gene name 
  gene_label <- data.frame(rownames( mat))
  colnames(gene_label) <- "primerid"
  rownames(gene_label) <- gene_label[, 1]
  #MAST requires a primerid (which essentially is the gene name)
  
   #create a MAST object, with the group id of each cell and the gene name 
  sca <- MAST::FromMatrix(
      exprsArray = mat,
      cData = group_label,
      fData = gene_label
  )

  #this is the differential expression model (called a hurdle model)
  #(see https://www.bioconductor.org/packages/release/bioc/vignettes/MAST/inst/doc/MAITAnalysis.html#4_differential_expression_using_a_hurdle_model)
  zlmCond <- MAST::zlm(~group, sca) 
  summaryCond <- summary(object = zlmCond, doLRT = 'groupgroup2') 
  
  #now extract the information 
  summaryDT <- data.frame(summaryCond$datatable)
  return_val <-data.frame(summaryDT[summaryDT[, "component"] == "H", 1], summaryDT[summaryDT[, "component"] == "H", 4])
  #We select the "H" rows, the "H" means we want hurdle P values
  #extract both the gene name and the P value associated with this gene
  colnames(return_val) <- c("gene", "P_value")
  
  #order by the p value , from the most significant 
  return_val <- return_val[order(return_val$P_value), ]
  
  return (return_val)

}

Here we provide an example of using the findmarker gene function.
To find out the marker gene in cluster 4, we type in 4 in the findmarker function. We then look at the top 10 genes ranked by p-value. We can then use ggplot to visualise the distribution of one of the genes across the dataset.

df<- data.frame(tsne_result$Y)
df<- cbind(df, as.factor(simlr_result_k6$y$cluster) )

df = df %>%
  dplyr::rename(
    tsne1 = X1,
    tsne2 = X2,
    cluster = `as.factor(simlr_result_k6$y$cluster)`
  )


marker <- findmarker(mat, simlr_result_k6, 4)

## Assuming data assay in position 1, with name et is log-transformed.

## 
## Done!

## Combining coefficients and standard errors

## Calculating log-fold changes

## Calculating likelihood ratio tests

## Refitting on reduced model...

## 
## Done!

marker[1:10,]

##          gene       P_value
## 12739    Gys2 6.068095e-153
## 19151    Shbg 2.513240e-134
## 12924     Hgd 1.577501e-126
## 19364 Slc27a2 8.229521e-124
## 6096     Gcgr 7.663681e-122
## 16209     Otc 1.338106e-115
## 12842      Hc 1.016114e-113
## 15057  Mogat2 1.509759e-113
## 12729    Gulo 6.801648e-113
## 14081   Lect2 1.976735e-112

index <- which(rownames(sce_scMerge[,ids]) == "Gys2")
# df_Gys2 <- cbind(df, mat[index,] )
df_Gys2 = df
df_Gys2$Gys2 = mat[index,]

ggplot(data =df_Gys2, mapping = aes(x=tsne1, y = tsne2, colour = Gys2) ) +
  geom_point(alpha=0.5) +
  scale_color_viridis() +
  theme_bw() +
  labs(col="Gys2 expression", x = "tsne1", y = "tsne2")

Quick quiz

1. What does the above t-SNE plot tell you?

Here we repeat the analysis as above, but for cluster 3. See if you can understand the output.

marker <- findmarker(mat, simlr_result_k6, 3)

## Assuming data assay in position 1, with name et is log-transformed.

## 
## Done!

## Combining coefficients and standard errors

## Calculating log-fold changes

## Calculating likelihood ratio tests

## Refitting on reduced model...

## 
## Done!

marker[1:10,]

##                gene       P_value
## 88    1700011H14Rik 1.187064e-210
## 5220         Erich5 8.248276e-200
## 1257       Adamts16 4.544699e-193
## 3539          Cldn3 2.663833e-177
## 20941        Tmem27 2.406824e-172
## 11554         Gm609 7.277614e-164
## 14560      Marveld3 3.901531e-152
## 21906          Vnn1 5.224090e-147
## 13356        Igfbp7 5.040702e-134
## 19960          Spp1 8.186524e-133

index <- which(rownames(sce_scMerge[,ids]) == "Erich5")
# df_erich<- cbind(df, mat[index,] )
df_erich <- df
df_erich$erich = mat[index,]

ggplot(data =df_erich, mapping = aes(x=tsne1, y = tsne2, colour = erich) ) +
  geom_point(alpha=0.5)+scale_color_viridis()+
  theme_bw()+
  labs(col="Erich5 expression")

Cell type proportions

plt.dat = data.frame(
  table(lab)
)

ggplot(plt.dat, aes(x = lab, y = Freq, fill = lab)) +
  geom_bar(stat = "identity") +
  labs(x = "Cell types", y = "Frequency", title = "Composition of cell types")

We observe that hepatoblast/hepatocyte is the largest population.